Meta-Analysis of Dose-Response Studies with the dosresmeta R package

Evidence Synthesis and Meta-Analysis in R, March 27-31, 2023

Alessio Crippa

Outline


A vignette of the dosresmeta package

Dose-response meta-analysis

Dose-response meta-analysis


Derive a dose-response curve from published results of multiple studies on the relation between a quantitative exposure (e.g. number of cigarettes, diet) and the occurrence of a health outcome (e.g. mortality risk, incidence of a disease).


Research questions:

  • Is there any association? What is the shape of the association?

  • What are the exposure values (doses) associated with the best or worst outcome?

  • What are the factors that can influence the dose–response shape?

Increasing number of citations


  • Several fields of application.

  • Many leading medical and epidemiological journals.

  • International health organizations and academic institutions.

  • Measures of public health impact.

Aggregated data

Motivating example

library(dosresmeta)
library(tidyverse)
data(coffee_mort)
gt(filter(coffee_mort, id < 3), groupname_col = c("id", "author"))
year type dose cases n logrr se gender area
1 - LeGrady et al.
1987 ci 0.5 57 249 0.0000000 0.0000000 M USA
1987 ci 2.5 136 655 -0.2876821 0.1391187 M USA
1987 ci 4.5 144 619 -0.1743534 0.1373198 M USA
1987 ci 6.5 115 387 0.0861777 0.1401409 M USA
2 - Rosengren et al.
1991 ci 0.0 17 192 0.0000000 0.0000000 M Europe
1991 ci 1.5 88 1121 -0.1203576 0.2531537 M Europe
1991 ci 3.5 155 2464 -0.3418342 0.2442559 M Europe
1991 ci 5.5 155 1986 -0.1261707 0.2440559 M Europe
1991 ci 7.5 39 575 -0.2665264 0.2784189 M Europe
1991 ci 9.5 24 427 -0.5108256 0.2802582 M Europe

Single study analysis



p_1s <- subset(coffee_mort, id == 1) %>% 
  mutate(
    low_rr = exp(logrr - qnorm(.975)*se),
    upp_rr = exp(logrr + qnorm(.975)*se)) %>% 
  ggplot(aes(dose, exp(logrr))) +
  geom_point() +
  geom_errorbar(aes(ymin = low_rr, ymax = upp_rr), 
                col = "grey", width = .1) +
  scale_y_continuous(trans = "log") +
  labs(x = "Coffee consumption, cups/day", 
       y = "Relative risk",
       title = "LeGrady et al., 1987")

  • Observations are not independent (same reference category).

  • The RR in the referent is 1 by definition (log RR = 0).

Trend analysis

Log-linear model: \(y = X \beta + \epsilon\)

\(y\) : vector of non-referent log \(RR\).
\(X\) : design matrix contains the assigned doses.

  • Model without intercept, so that \(RR = 1|x = 0\).
  • \(Cov(\epsilon)\) can be approximated depending on the study design (“cc” case-control, “ir” incidence-rate, “ci” cumulative-incidence).

Additional info: number of cases, and number of subjects or person-time.

covar.logrr(cases, n, logrr, se^2, type = "ci", 
            data = subset(coffee_mort, id == 1))
           [,1]       [,2]       [,3]
[1,] 0.01935402 0.01279718 0.01259433
[2,] 0.01279718 0.01885673 0.01254856
[3,] 0.01259433 0.01254856 0.01963947

Linear trend


Assuming a linear trend: \(y = \beta_1 x + \epsilon\)

mods_lin <- dosresmeta(formula = logrr ~ dose, se = se, type = "ci", 
           cases = cases, n = n, data = subset(coffee_mort, id == 1))
Epi::ci.exp(mods_lin)
     exp(Est.)      2.5%    97.5%
dose  1.033376 0.9914259 1.077101


Every 1 cup/day of coffee increase is associated with a 3% (95% CI: 0.99, 1.07) higher mortality risk.

Non-linear trend

library(rms)
k <- quantile(coffee_mort$dose, c(.1, .5, .9))
mods_spl <- dosresmeta(formula = logrr ~ rcs(dose, k), se = se, type = "ci", 
                       cases = cases, n = n, data = subset(coffee_mort, id == 1))
summary(mods_spl)
Call:  dosresmeta(formula = logrr ~ rcs(dose, k), type = "ci", cases = cases, 
    n = n, data = subset(coffee_mort, id == 1), se = se)

One-stage fixed-effects meta-analysis
Covariance approximation: Greenland & Longnecker

Chi2 model: X2 = 11.5464 (df = 2), p-value = 0.0031

Fixed-effects coefficients
                   Estimate  Std. Error        z  Pr(>|z|)  95%ci.lb  95%ci.ub    
rcs(dose, k)dose    -0.2048      0.0814  -2.5156    0.0119   -0.3644   -0.0452   *
rcs(dose, k)dose'    0.3930      0.1300   3.0225    0.0025    0.1382    0.6479  **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

1 study 3 values, 2 fixed and 0 random-effects parameters
 logLik      AIC      BIC  
 3.7764  -3.5527  -5.3555


  • The is evidence of a non-linear association (\(P\) = 0.0031, \(P\) non-linear = 0.0025).
  • The dose-response \(\beta\) coefficients are difficult to interpret.
  • Results can be presented in a graphical (or table) format.

Graphical presentation

xref <- .5
preds_spl <- data.frame(dose = c(xref, seq(0, 7, .1))) %>% 
  cbind(spl = predict(mods_spl, newdata = ., expo = T), lin = predict(mods_lin, newdata = ., expo = T))
p_1s + geom_line(data = preds_spl, aes(y = spl.pred, linetype = "Curve Spline")) +
  geom_line(data = preds_spl, aes(y = lin.pred, linetype = "Linear")) +
  labs(linetype = "Model")

Multiple studies analysis

p_all <- coffee_mort %>% mutate(coffee_mort, low_rr = exp(logrr - qnorm(.975)*se), upp_rr = exp(logrr + qnorm(.975)*se)) %>% 
  ggplot(aes(dose, exp(logrr))) + geom_point() +
  geom_errorbar(aes(ymin = low_rr, ymax = upp_rr),  col = "grey", width = .1) +
  facet_wrap(~ id, scales = "free") + scale_y_continuous(trans = "log", breaks = pretty_breaks()) +
  labs(x = "Coffee consumption, cups/day", y = "Relative risk")
p_all

Study-specific models

Apply the same dose-response analyses across studies.

modi <- split(coffee_mort, coffee_mort$id) %>% lapply(function(x) {
  list("lin" = dosresmeta(formula = logrr ~ dose, se = se, type = type, 
                          cases = cases, n = n, data = x),
       "spl" = dosresmeta(formula = logrr ~ rcs(dose, k), se = se, type = type, 
                          cases = cases, n = n, data = x))
})

And get model-based predictions.

predi <- lapply(unique(coffee_mort$id), function(i) {
  xref <- coffee_mort$dose[coffee_mort$id == i & coffee_mort$se == 0]
  data.frame(dose = c(xref, seq(0, max(coffee_mort$dose[coffee_mort$id == i]), .1))) %>% 
    cbind(spl = predict(modi[[as.character(i)]]$spl, newdata = ., expo = T), 
          lin = predict(modi[[as.character(i)]]$lin, newdata = ., expo = T))
}) %>% set_names(nm = unique(coffee_mort$id)) %>% 
  bind_rows(.id = "id")

p_all <- p_all + geom_line(data = predi, aes(y = spl.pred, linetype = "Curve Spline")) +
  geom_line(data = predi, aes(y = lin.pred, linetype = "Linear")) +
  labs(linetype = "Model")
p_all

Meta-analysis

Contrast and combine the dose-response coefficients \(b_i\).

\[b_i \sim N \left( \bar b , v_{b_i} + \tau^2 \right)\]

\[\hat {\bar b} = \frac{\sum_{i = 1}^K b_i w_i}{\sum_{i = 1}^K w_i}\] \(w_i = (v_{b_i} + \tau^2)^{-1}\)

  • \(\tau^2\) is the between-studies heterogeneity.

  • Test for overall association: \(H_0: \bar b = 0\).

  • Test (Cochran \(Q\)) and quantify the impact (\(I^2\)) of statistical heterogeneity.

Linear trend

bi_lin <- data.frame(
  id = names(modi),
  yi = sapply(modi, function(x) coef(x$lin)),
  vi = sapply(modi, function(x) vcov(x$lin))
)
head(bi_lin)
       id          yi           vi
1.dose  1  0.03283124 4.470869e-04
2.dose  2 -0.02360445 3.561441e-04
3.dose  3 -0.01430817 5.833618e-05
4.dose  4 -0.04777017 6.194806e-04
5.dose  5 -0.04736154 1.219069e-03
6.dose  6 -0.02027627 1.127937e-04


library(metafor)
rma_lin <- rma.uni(yi = yi, vi = vi, 
                   data = bi_lin)
summary(rma_lin)

Random-Effects Model (k = 22; tau^2 estimator: REML)

  logLik  deviance       AIC       BIC      AICc   
 37.5674  -75.1348  -71.1348  -69.0457  -70.4681   

tau^2 (estimated amount of total heterogeneity): 0.0003 (SE = 0.0002)
tau (square root of estimated tau^2 value):      0.0173
I^2 (total heterogeneity / total variability):   75.01%
H^2 (total variability / sampling variability):  4.00

Test for Heterogeneity:
Q(df = 21) = 77.0088, p-val < .0001

Model Results:

estimate      se     zval    pval    ci.lb    ci.ub      
 -0.0326  0.0051  -6.4424  <.0001  -0.0425  -0.0227  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Non-linear trend

bi_spl <- data.frame(
  id = names(modi),
  yi = unname(t(sapply(modi, function(x) 
    coef(x$spl)))),
  vi = t(sapply(modi, function(x) vcov(x$spl)))
)
head(bi_spl)
  id        yi.1       yi.2         vi.1          vi.2          vi.3        vi.4
1  1 -0.20483274 0.39304403 0.0066300437 -0.0102252526 -0.0102252526 0.016910322
2  2 -0.07391787 0.07231451 0.0071693533 -0.0097924956 -0.0097924956 0.014074567
3  3 -0.02578942 0.02177087 0.0005222641 -0.0008797053 -0.0008797053 0.001668107
4  4 -0.17681455 0.30128718 0.0053709253 -0.0110934653 -0.0110934653 0.025900538
5  5 -0.08529809 0.08507452 0.0107967876 -0.0214784900 -0.0214784900 0.048166538
6  6 -0.20065255 0.24137081 0.0034649927 -0.0044857506 -0.0044857506 0.006002614


library(mvmeta)
mvmeta_spl <- mvmeta(cbind(yi.1, yi.2 ), 
                     S = bi_spl[, c(4, 5, 7)], 
                     data = bi_spl)
summary(mvmeta_spl)
Call:  mvmeta(formula = cbind(yi.1, yi.2) ~ 1, S = bi_spl[, c(4, 5, 
    7)], data = bi_spl)

Multivariate random-effects meta-analysis
Dimension: 2
Estimation method: REML

Fixed-effects coefficients
      Estimate  Std. Error        z  Pr(>|z|)  95%ci.lb  95%ci.ub     
yi.1   -0.0811      0.0115  -7.0780    0.0000   -0.1035   -0.0586  ***
yi.2    0.1072      0.0217   4.9451    0.0000    0.0647    0.1497  ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Between-study random-effects (co)variance components
Structure: General positive-definite
      Std. Dev     Corr
yi.1    0.0372     yi.1
yi.2    0.0707  -0.9601

Multivariate Cochran Q-test for heterogeneity:
Q = 101.3912 (df = 42), p-value = 0.0000
I-square statistic = 58.6%

22 studies, 44 observations, 2 fixed and 3 random-effects parameters
  logLik       AIC       BIC  
 46.2962  -82.5924  -73.9040

The dosresmeta package

Why dosresmeta?


  • It facilitates both the study-specific dose-response and meta-analysis.

  • It gives a comprehensive framework for dose-response meta-analysis.

  • It provides functions for testing, predicting, and interpreting (presenting) results.

  • It implements extensions (one-stage procedure, continuous outcome).

Model estimation (linear)

mod_lin <- dosresmeta(formula = logrr ~ dose, id = id, se = se, type = type, 
                       cases = cases, n = n, data = coffee_mort, method = "mm")
summary(mod_lin)
Call:  dosresmeta(formula = logrr ~ dose, id = id, type = type, cases = cases, 
    n = n, data = coffee_mort, se = se, method = "mm")

Two-stage random-effects meta-analysis
Estimation method: 
Covariance approximation: Greenland & Longnecker

Chi2 model: X2 = 44.4226 (df = 1), p-value = 0.0000

Fixed-effects coefficients
             Estimate  Std. Error        z  Pr(>|z|)  95%ci.lb  95%ci.ub     
(Intercept)   -0.0324      0.0049  -6.6650    0.0000   -0.0419   -0.0229  ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Between-study random-effects (co)variance components
  Std. Dev
    0.0163

Univariate Cochran Q-test for residual heterogeneity:
Q = 77.0088 (df = 21), p-value = 0.0000
I-square statistic = 72.7%

22 studies, 22 values, 1 fixed and 1 random-effects parameters


Epi::ci.exp(mod_lin)
            exp(Est.)      2.5%    97.5%
(Intercept) 0.9681412 0.9589672 0.977403

Model estimation (non-linear)

mod_spl <- dosresmeta(formula = logrr ~ rcs(dose, k), id = id, se = se, type = type, 
                       cases = cases, n = n, data = coffee_mort, method = "mm")
summary(mod_spl)
Call:  dosresmeta(formula = logrr ~ rcs(dose, k), id = id, type = type, 
    cases = cases, n = n, data = coffee_mort, se = se, method = "mm")

Two-stage random-effects meta-analysis
Estimation method: 
Covariance approximation: Greenland & Longnecker

Chi2 model: X2 = 209.5069 (df = 2), p-value = 0.0000

Fixed-effects coefficients
                   Estimate  Std. Error        z  Pr(>|z|)  95%ci.lb  95%ci.ub     
rcs(dose, k)dose    -0.0764      0.0105  -7.3101    0.0000   -0.0969   -0.0559  ***
rcs(dose, k)dose'    0.1053      0.0229   4.5935    0.0000    0.0604    0.1502  ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Between-study random-effects (co)variance components
                   Std. Dev              Corr
rcs(dose, k)dose     0.0354  rcs(dose, k)dose
rcs(dose, k)dose'    0.0797                -1

Univariate Cochran Q-test for residual heterogeneity:
Q = 101.3912 (df = 42), p-value = 0.0000
I-square statistic = 58.6%

22 studies, 44 values, 2 fixed and 1 random-effects parameters

Graphical presentation

data.frame(dose = c(0, seq(0, 6, .1))) %>% 
  cbind(spl = predict(mod_spl, newdata = ., expo = T), lin = predict(mod_lin, newdata = ., expo = T)) %>% 
  ggplot(aes(dose, spl.pred)) + 
  geom_line(aes(linetype = "Curve Spline")) +
  geom_ribbon(aes(ymin = spl.ci.lb, ymax = spl.ci.ub), alpha = .1) +
  scale_y_continuous(trans = "log") +
  labs(x = "Coffee consumption, cups/day", y = "Relative risk", linetype = "Model")

Best Linear Unbiased Predictions (BLUP)

  • The normality assumption on the random-effects can be exploited to better predict study-specific curves.

  • The empirical Bayes estimates’ or ‘Best Linear Unbiased Predictions’ can be computed as:

\[\hat b_i = \hat {\bar b} + \tau^2 v_{b_i}^{-1}(b_i - \hat {\bar b})\]

  • BLUPs are a trade-off between study-specific and the average curve.
# study-specific dose-response coefficients (blup)
bi_blup <- t(apply(blup(mod_spl), 1, function(x) x + rbind(coef(mod_spl)))) %>% 
  cbind(id = as.double(levels(mod_spl$id)))
predi_blup <- lapply(unique(coffee_mort$id), function(i) {
  xref <- coffee_mort$dose[coffee_mort$id == i & coffee_mort$se == 0]
  newdata <- data.frame(dose = c(xref, seq(0, max(coffee_mort$dose[coffee_mort$id == i]), .1)))
  ctr.mat <- rcspline.eval(newdata$dose, k, inclx = T)
  ctr.mat <- scale(ctr.mat, center = ctr.mat[1, ], scale = FALSE) # important to center the variables at the ref
  data.frame(dose = newdata$dose, blup = exp(ctr.mat %*% bi_blup[bi_blup[, "id"] == i, 1:2]))
}) %>% set_names(nm = unique(coffee_mort$id)) %>% 
  bind_rows(.id = "id")
p_all + geom_line(data = predi_blup, aes(y = blup, linetype = "BLUP"), col = "blue")

Goodness of fit

  1. Deviance (D):
    • Total absolute distance between fitted and reported RRs.
    • Test for model specification.
  2. Coefficient of determination (\(R^2\)):
    • Descriptive measure of agreement.
    • Dimensionless measure bounded between 0 and 1.
  3. Plot of de-correlated residuals versus exposure.
    • Visual assessment of the goodness of fit.
    • Evaluate how the pooled dose–response curve fits the data by exposure levels.

gof(mod_spl)
Goodness-of-fit statistics:

Deviance test: 
D = 141.332 (df = 77), p-value = 0.000

Coefficient of determination R-squared: 0.679 
Adjusted R-squared: 0.671


filter(coffee_mort, se != 0) %>% bind_cols(gof(mod_spl)$tdata) %>% 
  ggplot(aes(dose, tresiduals)) +
  geom_point() + geom_smooth(se = FALSE) +
  labs(x = "Coffee consumption, cups/day", y = "Transformed residuals")

Meta-regression and stratified analysis

  • To further explain heterogeneity and/or explore differences in dose-response curves.

  • Prespecified study-level characteristics.

  • Limited number of data points: interpret with caution.

  • Stratified (or subgroups) analysis: divide the observations in groups and evaluate differences.

  • Meta-regression: formalization of stratified analysis, + investigation of continuous variables (be aware of pitfall).

Meta-regression

spl_regr <- dosresmeta(formula = logrr ~ rcs(dose, k), id = id, se = se, type = type, 
                       cases = cases, n = n, data = coffee_mort, mod = ~ area)
spl_regr
Call:  dosresmeta(formula = logrr ~ rcs(dose, k), id = id, type = type, 
    cases = cases, n = n, data = coffee_mort, mod = ~area, se = se)

Fixed-effects coefficients:
 rcs(dose, k)dose.(Intercept)  rcs(dose, k)dose'.(Intercept)     rcs(dose, k)dose.areaJapan    rcs(dose, k)dose'.areaJapan       rcs(dose, k)dose.areaUSA      rcs(dose, k)dose'.areaUSA  
                      -0.0898                         0.1014                        -0.0497                         0.1573                         0.0221                        -0.0073  

22 studies, 44 values, 6 fixed and 3 random-effects parameters
  logLik       AIC       BIC  
 45.9634  -73.9267  -59.1884


Is there any evidence of differential association across geographical areas?

waldtest(Sigma = vcov(spl_regr), b = coef(spl_regr), Terms = c(3:6))
Wald test:
----------

Chi-squared test:
X2 = 9.0, df = 4, P(> X2) = 0.06

expand.grid(dose = seq(0, 7, .2), area = levels(coffee_mort$area)) %>% 
  cbind(predict(spl_regr, newdata = ., expo = T)) %>% 
  ggplot(aes(x = dose, y = pred, col = area)) +
  geom_line() + 
  scale_y_continuous(trans = "log") +
  labs(x = "Coffee consumption, cups/day", y = "Relative risk", linetype = "Model")

Model extensions

One-stage dose–response meta-analysis for aggregated data

  • Unified approach which combines dose-response and pooling estimation (linear mixed models for summarized data).

  • Advantage:

    • Conceptually easier.
    • Fit more elaborate curves.
    • Avoid exclusion of studies with limited number of observations.
    • Allow a better comparison of alternative models.

Model fitting

data("coffee_mort_add")
coffee_mort_all <- bind_rows(coffee_mort, coffee_mort_add) %>% 
  arrange(id)


# reconstructing covariances
Slist <- split(coffee_mort_all, coffee_mort_all$id) %>% lapply(function(x) {
  if (any(is.na(x$cases) | is.na(x$n))) {
    diag(x$se[x$se != 0 & !is.na(x$se)]^2, nrow = sum(x$se != 0 & !is.na(x$se)))
  } else {
    covar.logrr(y = logrr, v = I(se^2), cases = cases, n = n, 
                type = type, data = x)
  }
})

# Splines without exclusion
spl_noexcl <- dosresmeta(logrr ~ rcs(dose, k), id = id, type = type,
                  cases = cases, n = n, se = se, data = coffee_mort_all, 
                  covariance = "user", Slist = Slist, proc = "1stage")

# Graphical comparison
data.frame(dose = c(0, seq(0, 6, .1))) %>% 
  cbind(spl_2s = predict(mod_spl, newdata = ., expo = T),
        spl_1s = predict(spl_noexcl, newdata = ., expo = T)) %>% 
  ggplot(aes(dose)) + 
  geom_line(aes(y = spl_2s.pred, linetype = "Two-stage")) +
  geom_line(aes(y = spl_1s.pred, linetype = "One-stage")) +
  scale_y_continuous(trans = "log") +
  labs(x = "Coffee consumption, cups/day", y = "Relative risk", 
       linetype = "Approach")

Dose-response meta-analysis of differences in means

Extension of the methodology where the outcome is expressed as means/differences in means instead of (log) relative risks.

References and codes

Codes


References

The R dosresmeta package